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Abstract 

We propose a new integer programming formulation for the problem of finding a maximum stable set 
of a graph based on representatives of stable sets. In addition, we investigate exact solutions provided by 
a Lagrangian decomposition of this formulation in which only one constraint is relaxed. Some computa- 
tional experiments were carried out with an effective multi-threaded implementation of our algorithm in 
• a multi-core system, and their results are presented. 
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> ' 1 Introduction 

With the emergence of multicore computers, methods of parallelizing performance-critical applications, par- 
ticularly combinatorial optimization problems, become a relevant issue. One major question is the identifi- 
cation of algorithmic methods amenable to parallel treatment in a shared memory setting. In this context, 
decomposition appears as a useful tool to benefit from the potential of the processing power made available 
' by many cores. The objective of this paper is to show a simple and effective decomposition method for a 

^) . classical combinatorial problem and its implementation in a multicore system. 
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1.1 The problem and its mathematical formulation 



Let G = (V,E) be an undirected, simple, nonempty and connected graph. We write n and m for \V\ and 
\E\, respectively. A stable set of G is a subset W C V such that, for every pair u and v of distinct vertices in 
W, uv £ E, where uv is the notation for the edge of G defined by u and v. Stable sets induce structures in 
graphs that model constraints over unrelated elements of a given set (represented by V) when some pairs of 
elements are in conflict (represented by E). In such situations, G can be seen as representing logical relations 
between elements such that every edge uv defines a valid inequality of the form 

x u + x v < 1, (1) 

being x u and x v binary variables associated with u and v, respectively. In other words, only vertices of a 
stable set of G are allowed to get 1 in a feasible solution [3]. 
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The canonical problem involving stable sets is related to the stability number of graphs. Given that 
W denotes the family of all stable sets of G, the stability number a(G) of the graph G is the maximum 
cardinality of a member of W [18]. The Stable Set Problem (SSP) consists of determining a(G), for a given 
(arbitrary) graph G. A number of structural results and algorithmic techniques exist in the literature for 
this problem (see, for example, [1] for its properties related to perfect graphs, [4, 11, 12, 15] for polyhedral 
properties, [3, 9] for its use for generating valid inequalities in general mixed integer problems, [22] for its use 
in a branch-and-cut algorithm, and [6, 19] for lift-and-project relaxations). The formulation usually adopted 
in these studies is 

a(G) = max{a;(y) | x e STAB(G)}, (2) 

where a; is a size-n binary vector indexed by the vertices of G, x(U), for U C V, stands for X^ugc/ x [ u \ an< ^ 
STAB(G) is the convex hull of the incidence vectors of stable sets (including the empty stable set) of G, also 
referred to as the stable set polytope of G. In spite of the fact that unstructured SSPs are difficult integer 
programming problems (indeed, it is NP-Hard to approximate a(G) within a factor n e for some e > [2]), 
the main aspect to be observed in connection with (2) is its ability to provide tight upper bounds when 
a partial description of STAB{G) is used. In these terms, the number of variables employed, in 0(n), is 
clearly an advantage. However, the number of constraints is a drawback that can be partially circumvented 
with the cutting-plane method for inequalities of the type 

x(U) < a(G[U]), where U C V and G\U] stands for the subgraph of G induced by U, 

referred to as rank inequalities [4, 22, 25]. 

1.2 Enumeration 

There are two enumeration approaches usually adopted to solve the SSP (or the equivalent problem of finding 
a maximum stable set in the complement of G, known as the Maximum Clique Problem) to optimality. In 
both of them, a search tree of subproblems is traversed. The first approach is the best-first branch-and-cut 
where the cutting-plane method is applied to determine upper bounds for the subproblems and the traversal 
of the search tree is conducted according to a non-increasing order of such bounds [22] . We concentrate our 
study in this approach motivated by two main facts. Firstly, best-first is well suitable for parallel treatment; 
secondly, since several other problems can be formulated as stable set problems, polyhedral and cutting 
techniques tend to have broader applicability [3, 9, 17]. 

The second enumeration approach is the combinatorial-based depth-first branch-and-bound obtained 
by introducing a pruning strategy in a backtrack enumeration. Computational experiments reported in the 
literature show that this approach can be efficiently implemented, with good performance when dense graphs 
are considered. To achieve this good performance in practice, the bounds used to guide the traversal of the 
tree are obtained with combinatorial greedy heuristics than run faster than polyhedral techniques [13, 21, 23]. 
The accuracy of such bounds highly depends on an ordering of the vertices. 

1.3 The reformulation and the decomposition 

The decomposition we propose for the SSP is based on a reformulation of the problem by means of the notion 
of representatives of stable sets. The use of representative vertices is an approach that appeared originally in 
the literature to formulate integer coloring problems [7, 8]. It has been used with other stable set problems 
successfully [10, 14]. It consists of partitioning the family of stable sets of G into n sub-families, each of 
them having a vertex of G as a representative. This simple technique allows the consideration of disjoint 
collections of stable sets almost independently. 

In this paper, we propose the use of representative vertices also to formulating the SSP in a way that 
leads to a natural decomposition of (2) into at most n subproblems, each one being an instance of the SSP 
associated with an induced subgraph of G. The SSP becomes the problem of finding the largest stable 
set represented by each vertex and, among these, taking the largest one. This decomposition potentially 
improves the performance of the cutting-plane method since it allows the application of separation heuristics 
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to graphs that are smaller than G, at the cost of increasing the number of variables by a factor that may be 
as large as n. 

The most important property of the reformulation is that the SSPs derived for each representative are 
loosely coupled, which makes it simple to attain our objective of decomposing the original problem into 
smaller independent subproblems. A Lagrangian relaxation is then described that incorporates a single con- 
straint in the objective function. This technique appeared originally in [16] for more unstructured problems. 
In the case studied in this paper, the decomposition of the relaxed problem becomes much simpler due to the 
structure of the representatives formulation. As a result, we are left with a number of smaller independent 
SSP subproblems, each of them associated with a certain subgraph of G. Some computational experiments 
are presented to show that this decomposition can be solved at reasonable computational cost for medium 
size graphs. In addition, results obtained with a multi-threaded implementation are reported, showing that 
the decomposition we propose is suitable to be naturally parallelized and executed in nowadays multi-core 
architectures. 

1.4 Organization of the paper 

The paper is organized as follows. The notation adopted is stated in Section 2. In Section 3, we present 
the representatives formulations of the SSP and an algorithm to solve the reduced subproblems generated 
by the Lagrangian decomposition. Finally, we close the paper with the experimental results (Section 4) and 
conclusions (Section 5). 

2 Preliminaries 

2.1 Notation 

The complementary graph of G is denoted G = {V, E), and the number of edges \E\ is represented by m. A 
clique of G is a stable set of G. A clique cover is a partition of V into cliques of G. In the definitions in the 
sequel, v is a member of V . A vertex u G V, u ^ V, is called a v's neighbor if uv G E. Otherwise, u and v 
are anti-neighbors. The notation N(v) stands for the neighborhood of v in G, i.e. the set {u G V \ uv G E}. 
On the other hand, the neighborhood of v in G, called the anti-neighborhood of v, is denoted by N(v). 

Let W be a subset of V. Write i£[W] for the set of edges of G[W], which in turn is the subgraph of G 
induced by W. In the particular case when W is N(v), then we use G(v) for G[N(v)]. 

An orientation of G is a mapping a : E — > V such that o~(uv) G {u, v}. Define the out-neighborhood of 
u as N + (u) = {v G N(u) \ o~(uv) = v} and its in-neighborhood to be N~(u) = {v G N(u) \ o~(uv) = u}. 
A vertex s G V is a source in a if o~(sv) = v, for all v G N(s). The orientation a is acyclic if every clique 
has exactly one source in a. If an orientation is given for G, the out- and in-anti-neighborhoods can be 
defined similarly and be denoted by N + (u) and N~(u), respectively. Additionally, we write G~(v) instead 
of G[N~(v)] and G+(v) for G[N+{v)}. 

2.2 Scaling STAB(G) 

The n-dimensional characteristic vector xw of a stable set W G W is such that xw[u] gets 1 if u G W, and 
if u G V \ W. Let W also denote the set constituted by the characteristic vectors of stable sets of G, i.e. 
by all vectors x G {0, 1}™ such that every edge inequality 

x[u] + x[v] < 1, uv G E, 

holds. STAB(G) is then given by the convex hull of vectors in W. Thus, each vector x G STAB(G) is a 
convex combination of the form 

x = V] Xwx w , X! Aw = 1, Xw G [0, 1], for all W G W. 

wew wew 
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The version of a set P scaled by p G K is obtained by multiplying the vectors of P by p (multiplying a 
vector by a scalar corresponds to multiplying every entry of the vector by this scalar). In this sense, the 
polytope STAB(p, G), obtained by scaling STAB(G) by p, is the convex hull of vectors in W scaled by p. 
Therefore, a vector x G STAB(p, G) can be expressed as 

x = ^2 ^wxw, ^2 = P, A w G [0,p], for all W G W. (3) 

It follows that every vector x G {0,p} n of STAB(p,G) satisfies the edge inequality, scaled by p, i.e. 

x[u] + x[v] < p, uv G -B. 

Special cases of scaled stable set polytopes are STAB(Q, G) = {0}, where denotes the null vector, and 
STAB{G) = STAB (I, G), which gives 

a(G) = Max x(V), subject to £ G STAB {I, G). 

3 Problem reformulation and its reduction 

Assume we are given an acyclic orientation a of G and a total order -< on V that respects er, which means 
that a(uv) = v if, and only if, u < v, for every uv G -B. Let the in- and out-anti-neighborhoods be defined 
according to a. 

3.1 Representatives formulation 

A way to describe a family {Wi, . . . , Wk} of k pairwise disjoint stable sets of G is to take r^, the minimal vertex 
of Wi according to -<, to be the representative of Wi that represents Wi\{ri} C A r+ (r i ), for alH € {1, . . . , k}. 
From this point of view, the family of stable sets is seen as a set of representatives {n, . . . , rk} and, for each 
representative n , a set of vertices (maybe empty) represented by n ■ Observe that Wi is the result of the 
union of its n with a stable set of G + (r;). This fact gives raise to a reformulation of the SSP, as follows. 
Define a binary variable x u , for all u G V, to indicate whether u is a representative, and a vector a;" of binary 
variables indexed by the vertices in N + (u) to indicate the vertices represented by u. With these variables, 
we write 

a(G)=Max ^ x u + x ll (N + (u)) (4) 

uev\T 

Subject to x u G STAB(x u , G + (u)), ueV\T (5) 

J2 x u < 1 (6) 

uev\T 

x„G{0,l}, ueV\T (7) 

The constraints (6)-(7) and the fact that STAB(0, G) = {0} assure that every non-null feasible solution x u , 
x u , for aU u G V \ T, contains exactly one representative, say r. In mathematical terms, ir r = 1 and x u = 0, 
for all u G V \ T, u ^ r. This leads all the scaled stable set polytopes in (5), except that of r, to be inactive, 
which means that they only contain the null vector. Since the scaled stable set polytopes are disjoint and 
the objective function is linear, the (non-null) optimal solution will be attained at a vertex of the unique 
active polytope STAB(x r , G + (r)), that is, at a characteristic vector of a stable set of G + (r). The objective 
function chooses r so that the stability number of G + (r) is maximum. 

The set T C V is composed by vertices that cannot represent a stable set larger than those represented 
by the vertices in V \ T. Such a set generalizes the notion of simplicial vertices used to fix variables in [20]. 
It can be obtained iteratively as follows. Given an edge uv in E[V \ T], the vertex v can be included in T 
if N + (y) C N + (u). This condition implies that if W is a stable set represented by v, then W — v + u is a 
stable set represented by u. 
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3.2 Lagrangian relaxation and its decomposition 

A simple inspection of the representatives formulation shows that it comprises disjoint SSP subproblems 
coupled by inequality (6). Thus, the SSPs become independent by simply bringing the coupling inequality 
in the objective function with a penalty A G R + . Defining this Lagrangian relaxation of the representatives 
formulation, we get the Lagrangian function 

L(x,X)= ]T (x u +x u (N+(u)))+X(l- x u)> 

u£V\T ueV\T 

for A > [5] . The resulting separable Lagrangian problem is given by 

a x = X+ a A' 

u£V\T 

where 

al = max{(l - X)x u + x u (N+ (u)) : x u e STAB(x u ,G + (u)),x u e {0, 1}} . (8) 

Note that = ma,x{a(G + (u)) + 1 — A, 0}, for every u € V \ T. This decomposition of the problem 
in subproblems, each of them being an instance of the SSP in a subgraph of G, can be solved with a 
combination of the Lagrangian dual descent method to compute A with the branch-and-cut method to solve 
each subproblem, as described next. 

3.3 Solving the dual Lagrangian problem 

In Lagrangian approaches, it is usual to define a dual problem, which consists of finding a configuration of 
the multipliers that optimize the bound obtained with the relaxation. For any A > 0, it is known that tt\ 
is an upper bound for a{G) [5]. In particular, if the Lagrangian multiplier A is a(G), then = 0, for all 
u E V\T, and therefore a\ = a(G). It follows that an optimum value of A is precisely a(G). Moreover, the 
maximum value of L(x, a(G)) is attained when x is such that X)«ey x u = 1- 

To compute the optimum multiplier, we use an iterative algorithm. Starting with an initial estimate A , 
our algorithm generates, at iteration £ > 0, a new estimate X l+1 depending upon the solution of a relaxation 
of a subproblem. The iterations are described in Algorithm 1, where the subproblem LP U mentioned at 
line 7 is a linear relaxation of (8) (we give more details of LP U in the next section). Unsolved subproblems 
are ordered in a non- increasing order of upper bounds in the priority queue Q. The branch-and-cut tree of 
each subproblem is traversed only once in a best-first order (lines 12-24) and the Lagrangian multiplier is 
updated simultaneously with the traversal of the trees. At iteration £, the second largest upper bound oi of 
the subproblems in Q gives the estimate A^ as indicated in line 13, which is used to guide the search in the 
sense that all nodes in a selected search tree whose values are greater than X e are explored (lines 17-24). 
This means that the branch-and-cut algorithm proceeds with a subproblem as long as the associated upper 
bound is large enough. 

An interesting characteristic of Algorithm 1 is that the estimate of the Lagrangian multiplier decreases 
monotonically until a large stable set is found at line 22. This starts a phase where the multiplier increases 
monotonically until its optimum value is attained. A variable x u is fixed at zero when the maximum stable 
set that can be represented by u e V is at most as large as the largest stable set found so far. 

Since the subproblems in Q are independent, line 18 is performed in parallel on different subproblems. 

4 Computational experiments 

Here we present the results of computational experiments carried out with the representatives formulation 
and the Lagrangian decomposition. In what follows, we briefly describe some details of the implementations 
and its comparison with the standard model based on the STAB{G) polytope (formulation (2)). One main 
interest of this comparison is the investigation of the tradeoff between the increase of the number of variables 
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Algorithm 1 Solving the Lagrangian relaxation 



Define a priority queue Q of subproblems, ordered according their upper bounds 
a <- 

for all u G V do 

a <— max{a, 1 + size of a maximal stable set of G + (u)} 
a u <— 1+ size of a clique cover of G + (u) 
if a u > a then 

Generate the linear program associated with G + (u) and call it subproblem LP U 
Add LP U to Q, associated with its current upper bound a u 



A <— second maximum upper bound of the subproblems generated above 
while Q ^ do 



A 



a' — a 

Q' <- {u e Q : A < a u } 

for all u G Q' do in parallel 

Remove subproblem LP U from Q 

while LP U is unsolved, a u > a, and A < a u do 

Select an unsolved node of LP U and apply the cutting-plane method to it 
Update a u according to the upper bounds of the remaining unsolved nodes of LP U 
if best known integral solution of LP U is greater than A then 
a <— value of the new solution 

Remove from Q all nodes with upper bound at most a 
if a u > a then 

Perform a branching of the branch- and- cut method associated with LP U 

for all u G Q 1 do in parallel 

if LP U is unsolved and a u > a then 

Add LP U to Q, associated with its current upper bound a u 

7T <- max{7r/2,0.05} 
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and the strengthening of the subproblems obtained by the decomposition. The second aspect investigated 
with the experiments is the performance of the parallel versions of the decomposition. 

In all implementations, the data structure used to store the graph is a 0-1 adjacency matrix. The 
implementation of formulation 2 consists of the construction of an initial linear relaxation and a standard 
branch-and-cut solver. To construct the initial linear relaxation, every edge of the graph generates an 
inequality of type (1). Moreover, a greedy algorithm is used to determine a clique cover of G. A rank 
inequality associated with each clique of such cover is added to the linear relaxation. 

Our implementations of (4)-(7) consists of the optimization of a core linear program (referred simply 
to as LP U in Algorithm 1) containing a partial description of STAB(G + (u)), for each u G V \ T. Further 
details of these implementations are given next. 

4.1 Pre-processing 

Before solving the formulations, our algorithm performs a pre-processing computation on G. One task 
accomplished by the pre-processing is the determination of an orientation of G, which is used to define the 
variables of the formulations. The orientation is based on an ordering of the vertices obtained with the 
following iterative algorithm. At each iteration I = 0, 1, . . . , n — 1, the £th vertex in the order is determined. 
Let V e denote the set of vertices that remain to be ordered after £ iterations. Naturally, V° = V. The 
vertex u to order at iteration I is chosen in such a way that the number of edges connecting u to vertices 
in V 1 is maximized. To break ties, the strategy proposed in [23] is used. In this strategy, a vertex whose 
anti- neighborhood induces a denser subgraph is ordered last. According to these criteria, the very first vertex 
of the ordering is one among those with maximum degree and, in addition, a maximal stable set is placed 
at the last positions of the order. In Table 1, some properties of the ordering obtained are shown. 

The second task of the pre-processing is to determine bounds for a(G), which is done based on the order 
on the vertices. We take as a lower bound the size of the maximal stable set defined by the last vertices in 
the order. This is accomplished by the ordering algorithm, without any increase in the time complexity The 
lower bounds so obtained are shown in the last column of Table 1. As an upper bound for a(G), a simple 
greedy heuristic, guided by the order on the vertices, is used to determine a clique cover of G + {u). This 
clique cover also gives upper bounds for the subproblems. For each u GV\T 7 we simply take the number of 
cliques intersecting N + (u) as the initial upper bound indicated at line 5 of Algorithm 1. The comparison of 
these upper bounds with the lower bound leads the pre-processing phase to discard those subproblems that 
do not contain any larger stable set. In Table 1, the number of valid vertices corresponds to the number of 
subproblems that are not discarded. There are two measures to evaluate the density of the graphs associated 
with the subproblems, namely: their maximum and average values. 

4.2 Handling the core linear programs 

The construction of the core linear program of a subproblem is also based on the order on the vertices, which 
gives the definition of the variables and initial constraints. The clique cover of G (that determining the 
upper bound for a(G)) is used to generate rank inequalities, added to the core linear program together with 
the corresponding edge inequalities. It should be noticed that this order is different from that one used to 
produce the clique cover used in the construction of the linear program associated with the implementation 
of the formulation derived from STAB(G). 

Each core linear program evolves from its initial state driven by a branch-and-cut algorithm, which 
has already been proved to be effective to this problem [17, 22]. The separation heuristic at line 18 in 
Algorithm 1 is implemented in the linear program solver used, CPLEX 11.0. The parameters (including the 
choice of separation heuristics to apply) used to obtain the results arc determined using the tunc facility of 
CPLEX 11.0 with the problem defined on the STAB(G). These are the parameters used for solving the 
Lagrangian decomposition. The times reported next do not consider the time required to initialize the core 
linear programs. 
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Instance 






Ordering properties 




G 


n 


a(G) 


Valid Max subgraph 


Av. subgraph 


Stable 




(density) 




vertices 


density 


density 


set size 






Random graphs 








g200.90 


200 (0.90) 


4.6 


168.0 


0.979 


0.855 


3.6 


g200 . 70 


200 (0.70) 


7.2 


186.4 


0.720 


0.616 


5.4 


g200.50 


200 (0.50) 


11 


187.8 


0.507 


0.407 


8.2 


g200 . 30 


200 (0.30) 


18 


183.6 


0.302 


0.221 


14.2 






DIMACS graphs 








c-fat200.1 


200 (0.92) 


12 









12 


c-fat200.2 


200 (0.84) 


24 


2 


0.253 


0.248 


24 


c-fat200.5 


200 (0.57) 


58 


7 


0.268 


0.249 


58 


c-fat500.1 


500 (0.93) 


14 


1 


0.265 


0.265 


14 


c-fat500.2 


500 (0.63) 


26 


1 


0.242 


0.242 


26 


c-fat500.5 


500 (0.96) 


64 


11 


0.243 


0.234 


64 


c-fat500.10 


500 (0.81) 


126 


28 


0.245 


0.234 


126 


brock200.2 


200 (0.50) 


12 


187 


0.512 


0.416 


7 


brock200.4 


200 (0.44) 


17 


186 


0.347 


0.258 


12 


brock400 . 2 


400 (0.25) 


29 


372 


0.255 


0.192 


19 


brock400.4 


400 (0.25) 


33 


379 


0.252 


0.192 


19 


p_hat300-l 


300 (0.76) 


8 


282 


0.720 


0.591 


8 


p_hat300-2 


300 (0.51) 


25 


273 


0.423 


0.241 


24 


p_hat300-3 


300 (0.26) 


36 


266 


0.243 


0.132 


31 


hamming8 . 4 


256 (0.36) 


16 


240 


1.000 


0.582 


1 


keller4 


171 (0.35) 


11 


154 


0.383 


0.387 


8 


san200 . . 7 . 2 


200 (0.30) 


18 


198 


1.000 


0.145 


1 


san200. 0.9.1 


200 (0.10) 


70 


186 


0.105 


0.035 


12 


san200 . . 9 . 2 


200 (0.10) 


60 


186 


0.099 


0.045 


12 


C125.9 


125 (0.10) 


34 


90 


0.101 


0.058 


33 


C250.9 


250 (0.10) 


[44,48] 


211 


0.101 


0.063 


37 


gen200 . 9 . 44 


200 (0.10) 


44 


166 


0.100 


0.057 


30 


gen200 .9.55 


200 (0.10) 


55 


159 


0.099 


0.059 


36 


mann . a27 


378 (0.01) 


126 


375 


1.000 


0.033 


1 


mann . a45 


1035(0.004) 


345 


1032 


1.000 


0.014 


1 



Table 1 : Properties of the graphs associated with the subproblcms 
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4.3 Results 



The algorithms have been implemented in Java and C on a dual quad-core based computer, with 4GBytes 
of memory and 3.0GHz of clock frequency. We used as test-bench some uniform random graphs and some 
instances extracted from DIMACS implementation challenge [24]. These instances are organized in Table 2 
according to the density of the graphs since it is the key parameter in the comparisons performed. The 
results reported are the following: the size of the graph (number of vertices n and density), the value of an 
optimum solution and the initial lower bound, and information about the two implemented approaches. For 
the implementation of formulation (2), the upper bound given by the fractional solution of the initial core 
linear program, the total number of branch-and-cut nodes explored, and the total running time (in seconds) 
for solving the core linear program. For the implementation of Algorithm 1, besides the total number of 
branch-and-cut nodes explored by all core linear programs and the corresponding running times, the initial 
value A of Lagrangian multiplier is also provided. The values shown for each random graph are the averages 
over 5 instances. The results for the parallel (multi-threaded) implementation were obtained using four 
concurrent threads. 





Instance 






STAB(G) 






Lagrangian 




G 




n 


a(G) 


LP 


Nodes 


Time 


A 


Nodes 


Time 




(density) 


(initial a) 






















Random graphs 










g200.90 


onn 


I n on\ 


4.6 (3.4) 




/IOC CiC\ 

428. OU 


9.92 


13.4 


400.20 


4.99 


g200.70 


onn 
zUU 


{ n vn\ 
[U. 1 U ) 


7.2 (5.4) 


1 CO 

21. OO 


oc a c en 


100.35 


24.6 


1251.2 


25.76 


g200.50 


200 


/'n cn\ 
(U.oO) 


11 (8.2) 


24.70 


15098.80 


208.38 


36.8 


11922.0 


64.02 


g200 . 30 


200 


(0.30) 


18 (14) 


35.65 


105922.4 


879.05 


52.4 


112019.4 


428.86 










DIMACS 


graphs 










c-f at200 . 1 


200 


(0.92) 


12 (2) 


15.00 


135 


2.81 








c-fat200.2 


200 


(0.84) 


24 (2) 


26.00 


35 


4.93 


28 


4 


1.92 


c-fat200.5 


200 


(0.57) 


58 (2) 


77.50 


21 


1.71 


59 


18 


0.92 


c-f at500. 1 


500 


(0.93) 


14 (2) 


17.00 


543 


68.08 


19 


251 


31.04 


c-fat500.2 


500 


(0.63) 


26 (2) 


30.91 


485 


175.32 


32 


254 


109.21 


c-fat500.5 


500 


(0.96) 


64 (2) 


74.75 


89 


178.01 


74 


118 


116.88 


c-fat500. 10 


500 


(0.81) 


126 (2) 


155.33 


13 


84.70 


127 


21 


31.68 


brock200.2 


200 


(0.50) 


12 (7) 


24.49 


4112 


105.8 


36 


7978 


93.29 


brock200.4 


200 


(0.44) 


17 (10) 


32.83 


36455 


393.4 


50 


38651 


200.41 


brock400.2 


400 


(0.25) 


29 (20) 


68.19 


317438 


> 18000 


102 


565318 


> 18000 


brock400.4 


400 


(0.25) 


33 (18) 


67.83 


261937 


> 18000 


102 


578937 


> 18000 


p_hat300-l 


300 


(0.76) 


8(7) 


27.00 


6181 


755.24 


29 


7382 


395.70 


p_hat300-2 


300 


(0.51) 


25 (13) 


45.13 


10288 


843.19 


57 


12664 


532.48 


p_hat300-3 


300 


(0.26) 


36 (30) 


56.65 


292443 


8681.0 


89 


286461 


3740.53 


hamming8 . 4 


256 


(0.36) 


16 (4) 


16.69 


1 


0.41 


33 


51 


2.23 


keller4 


171 


(0.35) 


11 (7) 


15.30 


5032 


28.73 


26 


1972 


12.37 


san200 . . 7 . 2 


200 


(0.30) 


18 (13) 


18.00 


1 


0.27 


40 


158 


1.52 


san200 .0.9.1 


200 


(0.10) 


70 (39) 


70.00 


1 


0.019 


89 


40 


0.59 


san200.0.9.2 


200 


(0.10) 


60 (28) 


60.00 


1 


0.033 


82 


40 


0.75 


C125.9 


125 


(0.10) 


34 (25) 


43.20 


3737 


5.97 


57 


4350 


2.78 


C250.9 


250 


(0.10) 


[44,48] (35) 


71.91 


2460234 


> 18000 


101 


4992455 


> 18000 


gen200 . 9 . 44 


200 


(0.10) 


44 (32) 


44.00 


1 


0.14 


74 


255 


2.21 


gen200 .9.55 


200 


(0.10) 


55 (27) 


55.62 


1 


0.072 


83 


40 


0.80 


mann. a27 


378 


(0.01) 


126 (117) 


135.0 


8272 


2.24 


241 


3191 


0.89 


maim. a45 


1035 


;o.oo4) 


345 (330) 


360.0 


56649 


84.24 


668 


124538 


118.85 



Table 2: Summary of results on some random and DIMACS benchmark graphs 
The results reflect the strength of the new formulation and the effectiveness of the Lagrangian decompo- 

9 



G 


Sequential 


4 cores 


2x4 


cores 




Nodes 


Time 


Work 


Speedup 


Work 


Speedup 






Random graphs 








g200.90 


400.20 


4.99 


1 


1.00 


1 


1.00 


g200.70 


1251.2 


25.76 


1 


1.01 


1 


1.01 


g200.50 


11922.0 


64.02 


1 


1.19 


1 


1.21 


g200.30 


112019.4 


428.86 


1 


1.96 


1 


2.31 






DIMACS graphs 








c-f at200. 1 










_ 


_ 


c-fat200.2 


4 


1.92 


1 


1.00 


1 


1.00 


c-f at200.5 


18 


0.92 


1 


1.02 


1 


1.05 


c-fat500. 1 


251 


31.04 


1 


1.00 


1 


1.00 


c-f at500.2 


254 


109.21 


1 


1.00 


1 


1.00 


c-f at500.5 


118 


116.88 


1 


0.99 


1 


0.98 


c-f at500. 10 


21 


31.68 


2.48 


1.00 


2.48 


1.02 


brock200.2 


7978 


93.29 


1 


1.61 


1 


1.62 


brock200.4 


38651 


200.41 


1 


2.65 


1 


2.92 


p_hat300-l 


7382 


395.70 


1 


1.12 


1 


1.02 


p_hat300-2 


12664 


532.48 


1 


1.00 


1 


1.00 


p_hat300-3 


286461 


3740.53 


1 


1.28 


1 


1.29 


hamming8 . 4 


51 


2.23 


1 


2.19 


1 


2.42 


keller4 


1972 


12.37 


1 


2.15 


1 


2.13 


san200.0.7.2 


158 


1.52 


1 


2.34 


1 


2.24 


san200. 0.9.1 


40 


0.59 


1 


1.48 


1 


2.19 


san200.0.9.2 


40 


0.75 


1 


2.08 


1 


2.50 


C125.9 


4350 


2.78 


0.99 


2.96 


1 


3.43 


gen200.9.44 


255 


2.21 


0.20 


3.03 


0.25 


3.11 


gen200.9.55 


40 


0.80 


1 


1.57 


1 


2.29 


mann. a27 


3191 


0.89 


1 


2.23 


1 


2.12 


mann. a45 


124538 


118.85 


1 


3.11 


1 


3.08 



Table 3: Parallel performance with 4 and 8 threads 
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sition. This new approach outperforms the implementation based on STAB(G) in almost all instances with 
density at least 30%, as long as the running time is considered. The only significant exception is hamming8 . 4, 
in which case several subproblcms giving the stability number as optimum solution arc solved once before 
being pruned. Naturally, the structure of the graph has great influence on the computation time. One of 
such relevant structural aspects is the number of maximum stable sets. In cases where this number is small, 
the decomposition given by the ordering of the vertices gives directly the optimum solution. For instance, 
for the graphs c-fat200.1 to c-fat500.10, the number of nodes explored in the branch-and-cut trees is 
smaller when the decomposition is used when compared to the formulation based on STAB(G). For low 
density graphs, our proposed decomposition is competitive in many cases, but pathological behaviors are 
observed in some DIMACS benchmark graphs. For instance, in the cases san200.0.7.2 to san200 . . 9 . 2, 
subproblems are almost as difficult as the original problem based on formulation (2). Another point that 
can be observed is an acceptable increase in the number of nodes explored, which is the case with several of 
the graphs c-f at, brock, and pjiat. 

The gain obtained with the multi-threaded parallel implementation is shown in Table 3, for 4 and 8 
threads. Two measures are used to evaluate the performance of this parallel implementation. First, the 
variation of the number of nodes of the branch-and-cut tree explored is given in the column entitled "work" . 
We observe no significant variation in this measure for all graphs, except for gen200 . 9 . 44, in which case the 
parallel execution explored much less nodes. Second, the speedup given by sequential time/parallel time is 
larger for denser graphs. In general, we observe that the Java mechanism used to schedule the threads does 
not incur any significant overhead. Another fact that should be noticed is that, for the graph mann. a45, the 
parallel implementation of the decomposition performs better than the implementation of ST AB(G). 

5 Concluding remarks 

A decomposition for a new formulation for the SSP is proposed and analyzed. An empirical comparison of 
this decomposition with the standard formulation shows a significant improvement in the efficiency of the 
branch-and-cut algorithm. According to experimental observations, the standard formulation performs well 
with sparse graphs. The decomposition leads to a better performance with denser graphs. Additional gain 
in performance is attained with a parallel implementation that explore the processing power of multicorc 
systems. 

As directions for further research, the effectiveness of the Lagrangian decomposition of the representatives 
formulation (and its shared memory parallel implementation) for the SSP makes this approach very promising 
for other problems involving stable sets, like the vertex coloring problem and its generalizations. 
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